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G\ • Abstract 

p i. We study the diffusion of helium and other heavy elements in the solar in- 
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terior by solving exactly the set of flow equations developed by Burgers for a 
multi-component fluid, including the residual heat-flow terms. No approximation 
1^ ! is made concerning the relative concentrations and no restriction is placed on the 



Oh' 



number of elements considered. We give improved diffusion velocities for hydro- 



O 

o 

^ ■ gen, helium, oxygen and iron, in the analytic form derived previously by Bahcall 

Q\ . and Loeb. These expressions for the diffusion velocities are simple to program in 



stellar evolution codes and are expected to be accurate to ~ 15%. We flnd that 



O ' the inclusion of the residual heat flow terms leads to an increase in the hydrogen 

in: 

. diffusion velocity. We compare our numerical results with those obtained ana- 

lytically by Bahcall and Loeb using a simplifled treatment, as well as with those 

•i-H . 

^ ' derived numerically by Michaud and Proffltt. We flnd that for conditions char- 

acteristic of the sun, the results of Bahcall and Loeb for the hydrogen diffusion 
velocity are smaller than our more accurate numerical results by ~ 30%, except 
very near the center where the error becomes larger. The Michaud and Proffitt 
results differ from the numerical results derived here by <15%. Our complete 
treatment of element diffusion can be directly incorporated in a standard stellar 
evolution code by means of an exportable subroutine, but, for convenience, we 
also give simple analytical flts to our numerical results. 

Subject headings: diffusion, stars: interiors, stars: abundances. Sun: int erior 
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1. Introduction 

Precise solar evolution calculations must be carried out to compare model results with 
observations of solar neutrino fliixes and of p-mode oscillation frequencies. In particular, 
element diffusion affects the element abundances, the mean molecular weight, and the 
radiative opacity in the core of the sun, and therefore affects the calculated neutrino fluxes 
and oscillation frequencies. The characteristic time for elements to diffuse a solar radius 
under solar conditions is of the order of 6 x 10^^ yrs, much larger than the age of the sun. 
Element diffusion therefore introduces only a small correction to standard solar model 
calculations. Bahcall and Pinsonneault (1992a,b) showed that helium diffusion increases 
the predicted event rates by about 11% in the chlorine solar neutrino experiment, by 3% 
in the gallium experiment, and by 12% in the Kamiokande and SNO experiments, while 
increasing the inferred primordial helium abundance by 0.4% and decreasing the calculated 
depth of the convection zone by 2%. Christensen-Dalsgaard, Proffitt and Thompson (1993) 
calculated the sound speed as a function of radius in the solar model and concluded that 
helium diffusion causes a significant difference in the computed radial profile of the sound 
speed. Guenther, Pinsonneault, and Bahcall (1993) demonstrated that helium diffusion 
has a characteristic which depends upon the degree and frequency of the p-mode being 
discussed and which has a typical amplitude of order 1-3 MHz. 

Since the effects of diffusion are small, there is in principle no need for very high 
accuracy in its treatment. However, discrepencies appear between various results in the 
literature, depending on the approximations made. Previous studies of element diffu- 
sion in the sun (AUer and Chapman 1960, Michaud et al. 1976, Noerdlinger 1977, 1978, 
Cox, Guzik and Kidman 1989, Paquette et al. 1986, Bahcall and Loeb 1990, Proffitt and 
Michaud 1991, Michaud and Proffitt 1992, Bahcall and Pinsonneault 1992a,b, Christensen- 
Dalsgaard, Proffitt and Thompson 1993, Guenther, Pinsonneault and Bahcall 1993, Vau- 
clair and Vauclair 1982 and references therein) have usually included one or more of the 
following simplifying assumptions: neglecting thermal diffusion, or treating it using a 
simplified empirical formula; neglecting the presence of heavy elements when calculating 
helium diffusion; assuming a negligible helium abundance when calculating the diffusion of 
heavier elements; adopting a single constant value for all Coulomb logarithms. In this pa- 
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per, we provide a simple but complete treatment of the problem, making none of the above 
approximations, and we compare our results with those obtained under different simpli- 
fying assumptions. In particular, we compare our results with those obtained by Bahcall 
and Loeb (1990) (hereafter BL) and those obtained by Michaud and Proffitt (1992) (here- 
after MP). BL made most of the above simplifying assumptions. In particular, they used 
empirical results for the thermal diffusion coefficients and a single value for the Coulomb 
logarithms, equal to 2.2. MP solved the Burgers equations and then represented the ef- 
fects of the residual heat flow vectors by an ad-hoc correction to the results obtained when 
neglecting those heat fluxes. The principal difference between this work and most previous 
studies is that we solve the Burgers equations exactly and then represent the numerical 
results by simple analytic functions, rather than trying to obtain analytic solutions by 
various approximations. 

Element diffusion in stars is driven by pressure gradients (or gravity), temperature 
gradients, composition gradients, and radiation pressure-*^. Gravity tends to concentrate 
the heavier elements towards the center of the star. In a pure hydrogen- helium plasma, 
helium diffuses towards the center of the star, while hydrogen diffuses outwards. As we 
will show in §4 (see also Bahcall and Loeb 1990), the local rate of change of the hydrogen 
mass fraction is equal and opposite to the rate of change of the helium mass fraction. 
This follows from the condition of momentum conservation. The light electrons also tend 
to rise, but are held back by an electric field which counteracts gravity. Temperature 
gradients lead to thermal diffusion, which tends to concentrate more highly charged and 
more massive species towards the hottest region of the star, its center. Concentration 
gradients oppose the above processes. Radiation pressure causes negligible diffusion in the 
solar core and will be neglected in this paper. 

We study the relative diffusion of hydrogen, helium, and heavier elements, such as 
oxygen and iron. In contrast to many previous studies, no approximation is made con- 
cerning the relative concentrations of the various species, and no restriction is placed on 

^ In this work, we ignore the effects of meridional circulation. It has been shown that 
meridional circulation velocities are several orders of magnitude smaller than the diffusion 
velocities in the solar interior (see, e.g., Michaud and Vauclair 1991). 



3 



the number of elements considered. Our method is therefore apphcable to a wide variety 
of astrophysical problems, such as the diffusion of elements in white dwarf envelopes (see, 
e.g., Fontaine and Michaud 1979, Pellcticr et al. 1986) and in globular cluster stars (see, 
e.g., Chaboyer et al. 1992). In this paper, we concentrate on calculating the diffusion ve- 
locities in the temperature and density ranges relevant to the sun, although our exportable 
subroutine can be used to calculate diffusion velocities in red giants and in white dwarfs. 

Burgers (1969) has provided a complete and straightforward set of equations to de- 
scribe the evolution of a multi-component fluid. In order to include the effects of thermal 
diffusion, he introduced the so-called "residual heat flow vectors". Here, we will use the 
Burgers equations, including the residual heat fluxes, to describe the plasma in the solar 
interior. Even though these equations can in principle be solved analytically, the algebraic 
complexity increases rapidly with the number of species considered. For example, because 
of computational limitations, Noerdlinger (1977) included only three species (hydrogen, 
helium, and electrons) and adopted a single constant for all the Coulomb logarithms. In 
contrast, we solve the full set of Burgers equations numerically, and place no restriction 
on the number of species. The Coulomb logarithm is obtained by calculating the collision 
integrals using a pure Coulomb potential with a long-range cutoff at the Debye length. 
However, the result obtained for the Coulomb logarithm is valid only for plasmas that are 
sufficiently hot and rarefied, i.e., such that the plasma parameter A is much larger than 
unity. For conditions characteristic of the solar interior, the Coulomb logarithms are small, 
and can even become negative for collisions between heavy elements. For such plasmas, the 
collision integrals can be calculated numerically using a screened Debye-Huckel Coulomb 
potential. The results can then be fitted to simple analytic functions. We adopt an ex- 
pression for the "effective" Coulomb logarithm obtained by Iben and MacDonald (1985) 
by fitting numerical results from Fontaine and Michaud (1979b). 

It should be relatively easy to incorporate our complete treatment of element diffusion 
into any standard solar evolution code^. However, we have obtained simple analytic fits to 
the exact results, which can provide a convenient alternative. These fits can be expressed 
as follows: Following BL's notations (see footnotes 3 and 5) and using BL's dimensionless 

^ Our FORTRAN routine will be made available upon request. 
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variables (see §2), the mass fraction of element s satisfies the equation 



ax. 



1 d 



(1) 



dt pr"^ dr 

where the partial derivatives are evaluated in the local rest frame of a mass shell in the 
star, i.e., in Lagrangian coordinates. The function ^s(r) is related to the diffusion velocity 
Ws of species s through 

Ur)^Ws{r)p{r)/T''/^{r). (2) 

We have obtained the following results for the diffusion velocities of hydrogen, oxygen and 
iron in the solar interior: 



, , , . , .dlnp ^ . ,ainT ^ , . SlnC^ 



dr 



dr 



dr 



(3) 



with 



(4) 



(5) 



' Ap{H) = -2.09 + 3.15 X - 1.07 X^, 

< ylT(iy) = -2.18 + 3.12 X- 0.96 X^ 

^ Ah{H) = -1.51 + 1.85X - 0.85X2, 
for the hydrogen diffusion coefficients, 

' Ap{0) = 0.15 + 1.34X - 0.89 X^ 

< At(0) = 0.53 + 1.99X- 0.72X2, 

[ Ah{0) = 0.08 + 0.58 X - 0.28 X^, 

for the oxygen diffusion coefficients, and 

Ap(Fe) = 0.25 + 1.31 X - 0.87X2, 

AriFe) = 0.65 + 1.99 X - 0.75 X^, 

^ AniFe) = 0.09 + 0.53X - 0.27X2, 

for the iron diffusion coefficients. These fits were obtained by using a constant value for 
each Coulomb logarithm, equal to its value at the center of the sun, and are accurate 
to better than 15% for the hydrogen and oxygen diffusion velocities, and better than 
20% for the iron diffusion velocity. Of course, the fits can only have a limited domain of 
applicability, whereas the numerical routine is completely general. 



(6) 



This paper is organized as follows. In §2, we introduce the notation and basic equa- 
tions. In §3, we describe the method of solution. In §4, we give the results for the hydrogen 
and helium diffusion coefficients, for a fixed value of the temperature and density, charac- 
teristic of the solar core. We compare these results with those obtained by BL and MP. 
In §5, we give the results for the heavy element diffusion coefficients, obtained under the 
same conditions. In §6, we give the diffusion velocities in the sun, and again we compare 
our results with those obtained by BL and MP. In §7, we give analytical expressions for 
our numerical results. In §8, we compare our expression for the electric field with the value 
obtained by Braginskii. Finally, in §9, we give a summary of the most important results. 

2. Basic Equations 

Each species of particles s is described by a distribution function Fs(x, v, t) normalized 
to unit integral, a mean number density Ug, an ionic charge Qs = ZgC, and a mass m^. 
All species are assumed to be at the same temperature T and in an overall hydrostatic 
equilibrium, since the temperature and pressure equilibration timescales are much shorter 
than the diffusion times. The mass and charge densities are ps = rigms and pes = '^s^s- 
The mean fluid velocity of each species is defined by 

Vis = j vFsdv. (7) 

The mean fluid velocity is given by 

u=-y'psUs, (8) 
P 

' s 

where p — ps is the total mass density. The diffusion velocity for species s is defined 

by 

Ws = Us - u, (9) 

and is therefore measured relative to the mean velocity of the fluid as a whole. We define 
the "residual heat fiow vector" for species s by (Burgers, 1969): 

(10) 



J Fs{v -u)\y -u\'^dv - 
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where A;^ is Boltzmann's constant. The cross-section for Coufomb scattering between 
particles of species s and of species t (s can be equal to t) is given by 

ast = 2V^e^Z^,Z^{kBT)-^lnAst, (11) 

where InA^^ is the Coulomb logarithm (A^^ is the "plasma parameter"), a correction factor 
taking into account the logarithmic contribution of binary encounters with impact param- 
eters up to the Debye shielding length. For the Coulomb logarithm, we adopt the following 
expression, obtained by Iben and MacDonald (1985) using numerical results from Fontaine 
and Michaud (197 9b), 



1 A 1-6249, 
InAst = — - — In 



1 + 0.18769 



4kBTX 



(12) 



where A = max(AD,ao), \d = j ^'Ke^^^ngZlYl'^ is the Debye length, and ao = 
(3/47r X^ions '^»)^^^ interionic di stance. The friction coefficient between species s 

and t is 

Kst = (2/3)/Xst(2/cBT//x,t)i/2 (13) 

where = msmt/{ms + rrit) is the reduced mass for species s and t. 

The Burgers equations for mass, momentum, and energy conservation can then be 
written as 

1 S(.V„.)=('^) . (14) 



nucl. 



dps 
dr 



dt dr \ dt 

+ psQ - pesE = ^ Kat {wt - Ws) + 0.6{Xstrs " yst^)] , (15) 



t^s 



and 

^"-s^B^ = ^Kst{^Xst{wa -wt)-yat l.Qxstijs + Tt) + Ystrs - 4..?>Xstr^ I - 0.8^, 

t^s 

(16) 

In these equations, g = — (GM(r)/r^)er is the gravitational acceleration, g = |g|, E is 
the electric field, E = |E|, Xst = l^st/ms, Vst = ^J-st/mt, and Yst = ^Ust + l.^Xstrnt/nis. 
The numerical coefficients in equations (15) and (16) ar e related to the collision integrals 
and were obtained using a pure Coulomb potential with a long-range cutoff at the Debye 
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length. More accurate results can be obtained by using numerical values derived from 
calculations using a sere ened Coulomb potential. We have assumed spherical symmetry 
and included a term for composition changes due to nuclear burning in equation (14). 
Using equation (15), it is straightforward to show that 

T.(^+Ps9-PesE]=0, (17) 

or 

^+pg-p^E = (18) 

where p = X^^Ps and pe = Pes ^ire the total pressure and total charge density. The 
departure from local charge neutrality is very small, with peE/pg ~ Grri^/e^ ~ 10~^^ 
(see discussion and eqs.(22)-(23) in BL). Equation (18) therefore reduces to the familiar 
equation of hydrostatic equilibrium, 

dp/dr = —pg. (19) 
In addition, the following constraints must be satisfied: charge neutrality, 

^ge^n, = 0; (20) 



current neutrality. 



and local mass conservation. 



qesrisWs = 0, (21) 
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^^msngWs = 0. (22) 



Note that equation (22) follows from the fact that the diffusion equations are solved in 
the rest frame of the plasma. The set of linear equations (15)- (16) and (21)- (22) forms a 
closed system for the diffusion velocities Wg, the residual heat flow vectors r^, the gravi- 
tational acceleration g and the electric field E in terms of the pressure, temperature, and 
concentration gradients. Since we already know the value of = GM(r)/r^, this relation 
provides a useful check on the numerical results. 



If we ignore thermal diffusion, the electric field is given by = —{l/en(.){dpf./dr). 
When thermal effects are included, the electric field can be written as (see eq. (57) in BL) 

eE = -^-aekBjr- 23 

Tie or or 

Using his two-component fluid equations, Braginskii (1965) has obtained the values ctg ~ 
0.71 for a pure hydrogen plasma and cte ~ 0.9 for a pure helium plasma. 

3. Method of solution 

The system of equations (15)-(16) and (21)-(22) can be solved numerically. If there are 
S species in the system {S — 1 ions plus electrons), there are S momentum equations (15) 
and S energy equations (16). The unknowns are the S drift velocities Ws and the S heat 
fluxes rs- The gravitational acceleration and the electric fleld are also treated as unknowns, 
and we use the two additional equations for mass and charge conservation, equations (21) 
and (22), to help determine g and E. Note that the hydrostatic equilibrium condition 
(eq. 19) should be satisfied automatically, providing a useful check on the numerical results. 

We now rewrite the basic equations in a dimensionless form that is better suited to a 
numerical treatment. The radius r is expressed in units of , the mass density p in units 
of 100 gcm""^ and temperature T in units of 10^ K, characteristic values at the center of 
the sun, and the time t is in units of tq = 6 x lO-'^'^yrs, a characteristic diffusion time in 
the sun (see,e.g., Kippenhahn and Weigert p. 60, or eq. (9) in BL). We write the Burgers 
equations (15)-(16) and the constraints (21)-(22) as 

S ,1 ^ 25+2 



P 



dlnp dliiT dlnCj ^ 



where the following notations have been introduced. The concentration of species s is 
defined by 

Cs = ris/rie. (25) 
It is related to the mass fractions Xg = rusris/p by 
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or inversely 



(27) 



where is the atomic number of species s, and the sum is over all species, including the 
electrons. For the electrons, 

Af, = m^lmQ, (28) 
where me and mo are the electron and atomic masses, and 



a = 1. 



(29) 



The constant Kq is given by 



Kq = 1.144 X IQ-^^ T-^/^nl, 



(30) 



where T and p are expressed in the characteristic units defined above. We use the ideal 
gas equation of state, Ps — UskBT, and equations (25), (28) and (29) to write 

J15/2 



^ = 2.00 



(31) 



where we have written the electron number density in terms of the mass density, p = 
TTiQUe Y^gAgCs- The variables Wi are 



Wi 



Ti-S 



for i = 1, ...S 
ioTi = S+1,...2S 



K^^rieeE for i = 2S + 1 



(32) 



Kq ^nemog for i = 2S + 2. 
If we define C = J2i ^^e coefficients on the left-hand-side of equation (24) are given by 



C 



for z = 1,2, ...5 



foTi = S + l,...2S + 2, 



„. _,2.5^^ for z = 5 + 1, ...25 

for z = 1,..., 5 and i = 25 + 1,25 + 2, 



and 



7ii 



9l 
C 




c 



C Z2C2 



for z = 1, 5 

for z = 5 + l,...,25+2. 



(33) 



(34) 



(35) 
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The coefficients on the right-hand-side of equation (24) are given by 



K 



for j = i 

for j — 1, S and j ^ i 
Y.k^i O-Q^ikXik ioT j = i + S 

-0.6Kij-syi,j-s for j = S + 1, 2S and j + S 



1-j 



(36) 



. —AiCi 



for j = 25 + 1 
for j = 25 + 2 



for i=l,...,S, by 



^k^j ^■5l^t-S,kXi-S,k 
-1.5Ki-S,jXi-S,j 

- Y^k^i '^i-S,kyi-S,ki^-QXi-S,k + Yi-S,k) 

-0.8Ki_s,i-s for j 

'^■7'^i~s,j-syi-s,j-sXi-s,j-s 
I 



for J = z — 5 

for j — 1, S and j ^ i — S 



(37) 



for j = S + 1, 2S and j 7^ z 
for j = 25 + 1, 25 + 2 



for i=S+l,...2S, by 



ZjCj for j = 1, 5 

for j = 5+ l,...,25 + 2 



(38) 



for i=2S+l, and finally by 



for j = 5 + 1, ...,25 + 2 



(39) 



for i=2S+2. In these expressions, the coefficient Kst is defined by 



AsAt 
As+At 



1/2 



C.CtZ^Z^lnA^t. 



(40) 



It is related to the friction coefficient through Kgt = KqKsi- We have used the constraint 
of charge neutrality to eliminate the concentration gradient of species 2 in equation (24), 



dlnC2 _ _ A ZjCj dlnCj 
-~ Z2C2 dr 



dr 



(41) 



Since p/Kq is proportional to T^/^/ p (see eq. 31), all the velocities will be proportional 
to T^l"^ I p. Therefore, we introduce the function (following BL), such that 



= (t^/Vp)6. 
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(42) 



The rate of change of the element mass fractions due to diffusion is now written in dimen- 
sionless form as 



^ ^ ^r'X,r^/%{r)], (43) 



dt pr'^ dr 

the generahzation of equation (1) of BL to the case of arbitrary concentrations and a more 
accurate treatment of the plasma physics. 

Equations (24) are linear. Therefore, wc can combine linearly the solutions obtained 
by keeping only one of the gradients different from zero. We write the results in terms of 
generalized diffusion coefficients Ap{s), At{s) and At{s) for species s, as 

. , , . , .dlnp . , . 91nT ^-^ ^ , ,91nCt 
Ur) = A,{s)^+At{s)^;^ + Ms)^^. (44) 

If InA is assumed identical for all the interactions, the coefficients Ap, At, and At are 
functions of the mass fractions only. If InA is defined by equation (12), these coefficients 
also depend on the charges, the temperature, and the density. 

4. Hydrogen and helium diffusion 

First, we consider the diffusion of hydrogen and helium, neglecting the presence of 
heavier elements. We calculate the hydrogen diffusion velocity. The helium diffusion 
velocity is then simply obtained from the constraint that there is no mean fluid velocity, 
J^s-^sWs = 0. Neglecting the electron mass compared to the proton mass, we have 
Wa = —{X/Y)wh- The rate of change of the helium number density is therefore equal 
and opposite to the rate of change of the hydrogen number density, (dY/dt) = —{dX/dt). 
Helium diffuses towards the center of the star, whereas hydrogen diffuses outwards. 

In the absence of heavy elements, the function is given by^ 

= MH)^ + At{H)^ + Ah{H)^—. (45) 



^ Note that BL define with the opposite sign. They also write in terms of the 
mass fraction gradient instead of the number concentration gradients. These are simply 
related by d\nCs/dr = dlnXs/dr-{^. ZiXi/Ai)-'^ '£.{ZjXj/Aj)dlnXj/dr or inversely, 
dlnXg/dr = dlnCg/dr — (^■AiCi)~^'^jAjCjdlnCj/dr. If hydrogen and helium are 
the only elements, dlnX/dr = (1 + X)d\n.CH / dr. 
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We have chosen hehum as element number 2, i.e., we write the diffusion velocity in terms 
of the hydrogen concentration gradient, using equation (41) to eliminate the helium con- 
centration gradient. 

Two major simplications are usually made when calculating the hydrogen and he- 
lium diffusion velocities in the absence of heavy elements. It is usually assumed that the 
Coulomb logarithms InAy are identical for all interactions. This allows the factorization 
of In A outside the function (see, eg., Noerdlinger 1977, Bahcall and Loeb 1990). In 
that case, the coefficients Ap, At and Ah depend only on the hydrogen (or helium) con- 
centration, not on density, temperature and ionic charges. The second simplification is 
to ignore the residual heat fluxes Vg. Then, the diffusion velocities are easier to calculate 
analytically, since there is no need for the heat equations (16) and the number of variables 
and equations is reduced by a factor of two. However, these simplifications can lead to 
large relative errors in the diffusion velocities. In particular it has been argued by MP that 
thermal diffusion can increase the diffusion velocities by 30%. 

In figure la-c, we show the variation of the coefficients Ap, At, and Ah with the 
hydrogen mass fraction X. To obtain these results, we have assumed T = 10^ K and 
p = 100 g cm~^, typical values in the core of the sun. The exact results are represented 
with solid lines, the results obtained neglecting the heat fluxes are represented with short- 
dashed lines, and the results obtained by keeping the heat fluxes but usin g InA = 2.2 
for all interactions^ are represented with long-dashed lines. If the heat fluxes are totally 
neglected, the two coefficients Ap and Ah are underestimated, and At = 0. In figure Id, 
we show the relative errors on the coefficients due to these approximations. The short- 
dashed line represents the error on Ap and Ah when the heat fiuxes are neglected. It 
can be as high as 45% for small values of the hydrogen mass fraction (not relevant to the 
sun). The long-dashed lines represent the errors due to In A = 2.2. The errors on Ap, Ah 
and At are smaller than 20%, except when X ~ 1. In the interior of the sun, X varies 
approximately between 0.3 and 0.7. For these values of X, the error does not exceed 20%. 



^ This value is usually considered representative of the Coulomb logarithms in the solar 
interior (see, e.g., Noerdlinger 1977 and BL). 
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4.1 Comparison with Bahcall and Loeb 

In order to keep the analytical calculations simple without neglecting the thermal 
effects one can use "effective" thermal diffusion coefficients (obtained through fits to the 
exact numerical results). In figure 2a, we show the ratio between the "exact" coefficients 
and those obtained by BL who neglected the residual heat fluxes, assumed a Coulomb 
logarithm of 2.2 for all the interactions, and used an "effective" thermal diffusion coefficient. 
The expressions obtained by BL are: 



v-RT/x .nrdlnp .pr^lnT .pr^lnCw 

^ (r) = <^ + ^ + ^ (46) 



with 



^ = -5(1 - X)/4, (47) 
A^^ = -6(1 -X)(X + 0.32)7(1.8 - 0.9X)(3 + 5X), (48) 



and 



A^^ = -{X + 3)/{3 + 5X). (49) 



The result for the thermal diffusion coefficient was obtained by fitting values obtained 
previously by AUer and Chapman (1960), Montmerle and Michaud (1976), and Noerdlinger 
(1978). For values of X between 0.3 and 0.7, the error made by BL on Ap and Ah is smaller 
than 40%, whereas the error on At is as large as 70%. However, as we will show in §6, 
large errors on At do not necessarily lead to large errors on the diffusion velocities, since 
the temperature gradient in the sun is much smaller than the pressure and concentration 
gradients. 

It is important to notice that the heat fliixes affect not only the thermal diffusion 
coefficients, but also the pressure and composition gradient coefficients. 
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and 



where 



4.2. Comparison with Michaud and Proffitt 

MP solved the Burgers equations with and without including the heat fluxes, then 
represented the effects of the heat fluxes by an ad-hoc correction to the results obtained 
when neglecting those heat fluxes. In our dimensionless variables, their results can be 
written as 

cMP(\ _ AMp9}nj) j^pdlnT MpdlnCH /^.^x 

with 

^ 4(0.7 + 0.3X)(lnA^j,/2.2)' 

.MP^_9 (1 (r.2\ 

^ 8 (0.7 + 0.3X) (In Aa.j//2.2) ' ^ ' 

aMP ^ {^ + ^) (ro\ 

^ (3 + 5X)(0.7+0.3X)(lnA^y/2.2)' ^ ^ 

In A,j, = -19.95 - ^ Inp - ^ In + ^ l^T. (54) 

In figure 2b, we show in solid lines the ratio of our coefficients and those obtained by 
MP. The difference between our results and those obtained by MP is smaller than 15% for 
At, and smaller than 5% for Ap and Ah- This small discrepency will be discussed in §6. 

5. Heavy element diffusion 

Because of the complexity caused by the addition of heavy elements, this problem has 
always been approached with additional simplifications (Vauclair et al. 1974, Noerdlinger 
1978, BL, MP). One common simplification is to assume a negligible helium concentration, 
therefore reducing the problem to a two-species situation. However, this assumption is not 
valid in the interior of the sun, where the characteristic mass fractions of hydrogen, helium, 
and oxygen are of the order of 0.34, 0.64 and 0.01 respectively (see, e.g., Bahcall 1990). 

The functions for the hydrogen and oxygen are now written as 

= Ap{H) +At{H) + Ah {H) + Ao{H) , (55) 
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and^ 



d Inp 



+ At{0) 



dr 



dlnCH 
dr 



+ ^o(0) 



a In Co 
dr 



(56) 
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It is interesting to show the variation of the diffusion coefficients as a function of the 
hydrogen mass fraction for fixed values of the temperature, density, and heavy element 
mass fraction. Indeed, as we will show in the next section, in the sun these parameters all 
vary with radius, and it is more difficult to extract the hydrogen mass fraction dependence 
itself. In figure 3, we show the four coefficients Ap[H), At{H), Ah{H) and Ao{H) as a 
function of the oxygen mass fraction X. These results were obtained using T = 10^ K, 
p = 100 g cm~^ and Z = 0.01, where Z is the oxygen mass fraction. These values are 
typical in the solar interior. Note that now X has a maximum value determined by 
Xmax = ^ — because of the charge neutrality constraint. We notice that the coefficient 
Ao{H) is two orders of magnitude smaller than the other coefficients. This was expected, 
since we have chosen Kho/ Kua Z/Y 10~^. The error made by neglecting the 
presence of oxygen when calculating the hydrogen diffusion velocity is smaller than 2%. 

In figure 4 we show the oxygen difi"usion coefficients Ap{0), At{0), Ah{0) and 
Ao{0). Again, the coefficient Ao{0) is much smaller than the other three coefficients, 
and can be neglected to the level of precision desired in these calculations. 



We now calculate the difi^usion velocities of hydrogen, helium, and heavier elements in 
the present solar interior (r < 0.7-Rq). We use values for the pressure, temperature, density 
and mass fractions of the contemporary sun, obtained from the standard solar model (table 
4.4, in Bahcall 1990). Since we have the radial profile of the pressure, temperature, and 
mass fractions, we can calculate their gradients. The coefficients Ap^ At-, Ah and Aq 
are computed using the tabulated values of T, p, X, y, and Z. The iron abundance is 
assumed to be uniform and given by \ogiQ{nFe/nH) = 6.82 — 12, and its ionization is 21. 

In figure 5, we show the radial variation of the hydrogen diffusion coefficients. In 
figure 6, we show the relative importance of the different terms in the hydrogen diffusion 



^ Note that BL define the function U such that dZ/dt = -(l/pr^) d[r'^XZT^/'^U/{'2- 
X)]dr. The function is related to the function through = Co(2 — X)/X. 



6. Diffusion velocities in the sun 
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velocity. The oxygen concentration gradient gives a negligible contribution to wh- The 
temperature term is not negligible, but is smaller than the pressure term (between 25% and 
50% of the pressure term). Therefore, a large error on the temperature diffusion coefficient 
does not necessarily lead to a large error for the diffusion velocities. 

The timescale for a change in the element abundances can be characterized by t = 
r/wH- To obtain the time t in units of to = 5 x lO^yrs, the age of the sun, we simply 
multiply the dimensionless time by to/ro, where tq is the characteristic diffusion time 
defined in § 3. The smaller the time t, the faster the element concentrations change. In 
figure 7, we show the variation of t with the radius. The fastest changes in the hydrogen 
concentration occurs at approximately 0.05i?Q, where tmin ~ 70 to- 

6.1. Comparison with Bahcall and Loeb 

In figure 8a, we compare our exact results with those obtained using the analytic BL 
formulae, equations (47)- (49) (eqs. 1-5 in BL). The BL formula underestimate the diffusion 
coefficients. The error on the pressure and concentration diffusion coefficients is smaller 
than 30%. The error on the temperature diffusion coefficient is of the order of 50%, except 
near the center where it becomes very large. 

In figure 9a-c, we show in solid lines the results for the diffusion velocities of hydrogen, 
oxygen, and iron. The helium diffusion velocity is related to the hydrogen and oxygen 
diffusion velocities through the zero mean velocity constraint, equation (22), 

Wa = -{XhWh + XqWo + XFeWFe)/XHe- (57) 

We also show in short-dashed lines the BL results for the hydrogen and oxygen diffusion 
velocities, given by equations (47)-(49) for the hydrogen velocity, and equations (2)-(5) 
in BL for the oxygen velocity. In the BL approximation, the helium velocity is given by 
Wa = —{Xh / XHe)wH- The error in wh due to the BL approximations is smaller than 
30%, except near the center, where the error is as high as 70%. However, the error for the 
oxygen diffusion velocity is very large. This was expected, since BL neglected completely 
the presence of helium when calculating wq- The BL results for wq are therefore only 
valid when X 
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6.2. Comparison with Michaud and Proffitt 



In figure 8b, we compare our results with those obtained using the MP formula, 
equations (51)-(53). The difll^erence between our results and the results obtained by MP 
for the diffusion coefficients is smaller than 5% for Ap and Ah- Our thermal diffusion 
coefficient At is larger by about 20%. This discrepency may result from the fact that we 
have used Zst = 0.6 for the heat fiux coefficien t in equations (15) and (16). If the collision 
integrals are calculated using a screened Coulomb potential, the value obtained for Zst is 
about 2/3 smaller (Proffitt 1993). 

In figure 9a-c, we show the results for the diffusion velocities of hydrogen, oxygen and 
iron. The MP result for the heavy elements difi'usion velocities are given by 



Xj ( i+x V 

5X+3 \^5X+3 ) 



+ [l + Z,-A,{^)]^ 



-^i^ix ^ix ^iy Ciy) + A^y 



X{A]i'C,,-A](;c,y) + A]l''c,y 



-0.23 



(58) 



+ 



0.54(4.75X + 2.25)dlnT 



(In Aa;j, + 5) dr 

where Aij = AiAj/ [A^ + Aj) is the reduced mass in atomic number units, 
Cij = ln[exp(1.21nAy) + 1]/1.2, and InA^ is given by equation (54). The difference 
between our results and those obtained by MP is smaller than 15% for the hydrogen 
diffusion velocity, and smaller than 20% for the oxygen diffusion velocity. 



7. Analytical fits 

All the results shown above were obtained using the expression (12) for the "effective" 
Coulomb logarithms. If we use (as in the previous results) the correct charge, temperature 
and density dependent expression for the Coulomb logarithms, we cannot give a simple 
analytical fit, as in BL, for the diffusion coefficients in terms of the hydrogen mass fraction. 
Even though one could in principle incorporate a subroutine which solves the problem of 
element diffusion in a standard solar model evolution code, it is useful to give a simple 
analytical fit of the results obtained here in the solar interior. It is convenient to provide 
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to stellar evolution programs other types of input physics, such as opacities or equations 
of state, in the form of tabulated values, or in terms of approximate analytical fits. In 
order to provide a similarly convenient service for diffusion, we adopted a constant for 
each Coulomb logarithm, the value it has at the center of the sun. In figure 10, we show 
the values of the various Coulomb logarithms in the sun computed using equation (12). In 
figure lla-c, we compare the diffusion coefficients obtained with these Coulomb logarithms 
with those obtained using a constant value of A (equal to its central value). In figure lid, 
we show that the error on wh is smaller than 4% in the solar core {r/R^ < 0.4), and 
remains smaller than 15% up to the convection zone. The error on the heavy elements 
diffusion velocities are ^6% in the solar core, and remain <20% up t o the convection 
zone. We can now fit these results to second order polynomials. Since the presence of 
oxygen and iron have very little influence on the hydrogen velocity, we can assume that 
the diffusion coefficients depend only on X. We find: 

' Ap{H) = -2.09 + 3.15X - 1.07X^ 

< AT(ii') = -2.18 + 3.12 X- 0.96X2, ^gg^ 
^ Ah{H) = -1.51 + 1.85X - 0.85X^ 
for the hydrogen diffusion coefficients, 

' Ap{0) = 0.15 + 1.34X - 0.89X2, 

< ^t(0) = 0.53+ 1.99X- 0.72X2, ^gQ^ 
^Ah{0)^ 0.08 + 0.58 X - 0.28 X^, 
for the oxygen diffusion coefficients, and 

Ap{Fe) = 0.25 + 1.31 X - 0.87X2, 

AriFe) = 0.65 + 1.99 X - 0.75 X2, (61) 

AniFe) = 0.09 + 0.53X - 0.27X2, 

for the iron diffusion coefficients. The errors due to the polynomial fits are of the order 
of 0.2%. There is no need to use higher order polynomials, since the error made with 
these second order polynomials are already much smaller than the errors introduced by 
the simplified Coulomb logarithm. 
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For an easy comparison with BL and MP results, we can factorize the BL results and 
write the numerical fits as 



A,{H) = --{1-X)a„ 



At{H) 



and 



Ah{H) = 



6(1 + 0.32) 
(1.8-0.9X)(3 + 5X) 

(X + 3) 



at, 



-an- 



(62) 
(63) 

(64) 



(3 + 5X) 

First order polynomial fits give the following analytical results: 

r ap = 1.66 - 0.82 X, 

OT = 4.46- 3.65 X, (65) 
^ an = 1.63- 0.74 X. 

The error introduced by the polynomial fit is much smaller than the error due to the sim- 
plified Coulomb logarithm. Equations (62)- (65) can be used to improve existing diffusion 
subroutines that are based on the BL formalism. 

It is important to remember that these fits are made for the standard parameters of 
the solar interior. 

8. Hydrostatic equilibrium and electric field 

As explained in § 2, we solve the system of equations (24) for the diffusion velocities, 
the residual heat flow vectors, the gravitational acceleration, and the electric field. On 
the other hand, the equation of hydrostatic equilibrium given by equation (19) must be 
satisfied. The numerical results for g can be written as 



9 = 



, , , (i In » . , .d\nT , , .din Ch 



dr 



dr 

and Ar 



dr 



(66) 



In hydrostatic equilibrium, At = Ah = and Ap = 1. The equation of hydrostatic 
equilibrium is satisfied to 10~^ by the numerical results. 
Similarly, we write the results for the electric field as 



E 



-Pe 



eup 



Ap{E)—^+At{E)—j—+Ah{E) 



dr 



dr 



dr 



(67) 
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The coefficients in equation (67) should be compared with equation (23). We find an 
excellent agreement between our numerical results and equation (23) for the coefficients 
in front of the pressure and concentration gradients (with an error smaller than 0.2 x 
10~^). Our values for the coefficient At{E) are slightly larger than the results obtained 
by Braginskii (1965) for ae- As shown in figure 12, we observe however the same tendency 
for this coefficient to increase slightly with an increase in the helium concentration. We 
get CKe = 0.8 for X = 1, and ae = 0.9 for X = 0. 

9. Summary and conclusions 

We have developed a Fortran program to solve numerically the Burgers equations for 
an arbitrary number of species, without any approximation. For the discussion of solar 
conditions given here, we have neglected the radiative forces, but these forces could easily 
be incorporated in the numerical routine. The accuracy of the results for the diffusion 
velocities is limited only by the validity of the expression for the Coulomb logarithm. The 
diffusion velocities of hydrogen, oxygen and iron were calculated in the solar interior and 
compared with the results of Bahcall and Loeb (1990) and by Michaud and Proffitt (1992). 
The results of BL for the hydrogen diffusion velocity are smaller by ~ 30%, except near 
the center, where the error is much larger. The results obtained by MP for the hydrogen 
and oxygen diffusion velocities differ by ^15%. 

We provide analytical fits of our numerical results for the diffusion coefficients as a 
function of the hydrogen mass fraction in the solar interior (eqs. 59- blf). These fits were 
obtained by assuming fixed values for the Coulomb logarithms, equal to their values at the 
center of the sun. 
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FIGURE CAPTIONS 



FIG. 1. — Variation of the hydrogen diffusion coefficients with the hydrogen mass frac- 
tion in a pure hydrogen- hehum plasma, with T = 10^ K and p ~ lOOgcm"^. The soHd 
hues represent the results obtained using Burgers equations, with no approximation. The 
short-dashed lines represent the results obtained when neglecting the heat fluxes. The long- 
dashed lines are the results obtained by using a single value for all the Coulomb logarithms, 
equal to 2.2. The dash-dot lines are the results obtained when neglecting the heat fluxes 
and using a single value for all the Coulomb logarithms, equal to 2.2. (a) Pressure gra- 
dient coefficient Ap. (b) Temperature gradient coefficient At- (c) Hydrogen concentration 
gradient coefficient Ah- (d) Relative errors due to the approximations. The short-dashed 
line is the error on Ap and Ah when the heat fluxes are neglected. The long-dashed lines 
are the errors made by using InA = 2.2, and the dash-dot line is the error on Ap and Ah 
when both approximations are made. 

FIG. 2. — (a) Ratio of the exact hydrogen diffusion coefficients and those obtained by 
various approximations, in terms of the hydrogen mass fraction X, with T = 10^ K and 
p = lOOgcm""^. (a) Comparison with Bahcall and Loeb (1990). (b) Comparison with 
Michaud and Proffitt (1992). The solid lines and the dashed lines are the results obtained 
using equations (54) and equation (12) respectively for the Coulomb logarithms in the 
Michaud and Proffitt formulae (51)-(53). 

FIG. 3. — Variation of the hydrogen diffusion coefficients with the hydrogen mass fraction 
in a hydrogen-helium-oxygen plasma, with T = 10^ K, p = lOOgcm"^, and Z = 0.01. 

FIG. 4. — Variation of the oxygen diffusion coefficients with the hydrogen mass fraction in 
a hydrogen-helium-oxygen plasma, with T = 10^ K, p — lOOgcm"^, and Z = 0.01. 

FIG. 5 — Hydrogen diffusion coefficients in the present sun, as a function of radius. 
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FIG. 6 — Contributions to the hydrogen diffusion velocity in the sun due to each gradient. 

FIG. 7 — Local diffusion time of hydrogen, oxygen, and iron as a function of radius, in 
units of the age of the sun. 

FIG. 8 — (a) Ratio of the exact hydrogen diffusion coefficients and those obtained by 
Bahcall and Loeb (1990) as a function of the radius, (b) Ratio of the exact hydrogen 
diffusion coefficients and those obtained by Michaud and Proffitt (1992) as a function of 
the radius. 

FIG 9 — Diffusion velocities in the contemporary sun. The solid lines are the exact results; 
The short-dashed lines are the results of Bahcall and Loeb (1990) (see eqs. 1-5 in BL); The 
long-dashed lines are the results of Michaud and Proffitt (1992). (a) Hydrogen diffusion 
velocity, (b) Oxygen diffusion velocity, (c) Iron diffusion velocity. 

FIG. 10 — Coulomb logarithms in the sun. 

FIG. 11 — Diffusion coefficients in the sun, as a function of radius. The solid lines rep- 
resent the results obtained using Burgers equations, with no approximation. The dashed 
lines represent the results obtained using InA = 3.2 for all the interactions, (a) Pressure 
gradient coefficient Ap. (b) Temperature gradient coefficient At- (c) Hydrogen concen- 
tration gradient coefficient Ah- (d) Relative error on the hydrogen and oxygen diffusion 
velocities made by using InA = 3.2. 

FIG. 12 — Thermal coefficient for the electric field (see eq. 23) 
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